This document is a template for analysis and reporting eDNA results. This file is a Quarto document (a newer version of a Rmarkdown file Intro to a Quarto doc) that can be opened in Rstudio and serve two purposes: 1) run code and 2) build a html (or pdf) report document.
The authors intended this file along with the other folders & files within this project to document analysis steps, make analyses reproducible, and make future analyses smoother and quicker to run. See the Quarto Website for a more in depth introduction to Quarto docs. The “Document Purpose & Description” section can be deleted when creating a new analysis.
1.1 Document Components
YAML (yet another markdown language) header - the code contained at the top of this document bounded by a set of --- marks. This is a place to store document metadata and instructions for how the rest of the quarto doc should run code chunks and format rendered documents. See here for more information –> Authoring Docs and YAML HTML Options. We currently have the YAML set so when you Render this document a self-contained .html document will be created, this file format and formatting options can be further customized.
Code chunks - code chunks are like mini code scripts. In a quarto doc you can run these individually by clicking on the green arrow/triangle in the upper right corner of the chunk OR using your “run code” shortcut. These chunks are allow to user to write text in between code, and provide instructions on how the code should be incorporated into the rendered document, including how figures should be displayed. See here for more information –> Code chunk options (work in both quarto and Rmarkdown docs)
Markdown - This bullet list is written in markdown. Markdown is a simple text and formatting language. Use this to write text in the rendered report document. Here is a quick guide to markdown syntax –> Markdown cheat sheet. You can also use values and/or objects generated in your code to populate in your markdown text, see here —> How to use a code value in markdown text
1.2 Code Chunk Example - Run Code
The primary function of this document is to run code. Below is an example of a code chunk that runs R code. These are where data cleaning and analysis steps are contained. In Rstudio you can run code several ways:
Code
# Code chunk example# The coding language is indicated in the very first line, within the curly brackets {}. # Use `#| ` to set options for the individual code chunk below the first line of the code chunk.# A code chunk execute option `echo: true` enables the printing of code output in rendered report # here is example R code. To write comments within a code chunk, use a # at the beginning of the line, lines with # will not be run within a code chunk.test <-1+1test_2 <- test *6test_2
[1] 12
line-by-line (just like a normal .R file)
chunk-by-chunk (press the green right-facing triangle at the top right corner of the code chunk)
all code top to bottom ()
1.3 Render Quarto Docs - Build Report Doc
The secondary function of this document is to build a report document. This requires all the code in the Quarto document to run cleanly without errors. To build the report, click the Render button shown below. This will start a multistage process that runs all the code and markdown sequentially and formats it based on the YAML header and code chunk execute options.
There are two ways to alter/customize the resulting report:
The YAML header (described above). For example, to set the default execution options for all code chunks, you can add options in the YAML (execute:)
The final template .html product will look something like this:
2 Analysis Template
Below is the code required for eDNA analysis
2.1 Setup
Load packages
We used the renv package to manage & record package versions. This records the version of R, package versions, and their dependencies into a file renv.lock at the root of this .Rproj. The purpose of this is to allow users on different computer setups and in the future to reproduce the results of this project in the same way. This is because different package versions can behave differently version to version, and even introduce errors or produced different results. By using the same exact set of package versions that the authors used to write the code, the same results can be reproduced and prevent unintended errors.
Renv creates a separate custom project library that is separate from the system library typically used when the library() function is called. This means that when you recreate the renv project library by calling renv::restore() you will likely need to install many packages. Many may be already installed in your system library, but they will need to be installed again into this custom renv project library with the exact version recorded in the renv.lock file. All files, folders, and scripts within the .Rproj will share the renv project library.
For more information about renv and how to add/update the renv project library -> .
Warning
Although this package management strategy is meant to reduce errors among future users, it is not perfect. The authors found that the renv project library could not yet be recreated on a M1 chip Apple computer (2023-10).
Create paths to folders
Set analysis date
Load data
2.2 Clean & save data
Clean & organize data
The code chunk below “clean-data” is separated into its own ./scripts/clean_data.R file and sourced by calling source(file.path(dir_scripts, clean_data.R) here. If you are looking save space in this analysis document, or want to reuse this exact same data cleaning protocol in different analysis document. The same goes for “create-nearshore-offshore-data” and “create-near-offshore-species-list” chunks, they could be their own scripts or be added to a clean_data script. If you choose to do this, the “example-source-clean-data” chunk below can be used.
loaded exsisting clean data
Create Nearshore & Offshore dataframes (or load if they already exist)
Create Nearshore & Offshore species dataframes (or load if they already exist)
2.3 Data Analysis
Unique Species
Unique Species - nearshore
Unique Species - offshore
Unique Genera
Unique Families
Transform data to Phyloseq Object
Explore Components of Phyloseq Object (if needed)
Subset Phyloseq Samples - Keep Environmental Samples Only
Subset Phyloseq Environmental Samples - ESVs with marine fish and mammal annotations
Environmental Marine Samples With More Than 10,000 Sequences
Standardize number of reads to the median sequencing depth
---title: "Palau eDNA Report Template"subtitle: "Subtitle Here"author: - name: "Author name 1" affiliation: "Organization 1" email: "author1@email.com" - name: "Author name 2" afflication: "Organizaton 2" email: "author2@email.com"date: "`r Sys.Date()`"format: html: number-sections: true toc: true code-tools: true theme: mintly number-depth: 2 self-contained: truetitle-block-banner: "#33ccffff"title-block-banner-color: "#ffdd00ff"code-fold: trueeditor: visualeditor_options: chunk_output_type: inlineexecute: eval: true echo: false warning: false---# Document PurposeThis document is a template for analysis and reporting eDNA results. This file is a Quarto document (a newer version of a Rmarkdown file [Intro to a Quarto doc](https://quarto.org/docs/faq/rmarkdown.html){.uri}) that can be opened in Rstudio and serve two purposes: 1) run code and 2) build a html (or pdf) report document.The authors intended this file along with the other folders & files within this project to document analysis steps, make analyses reproducible, and make future analyses smoother and quicker to run. See the [Quarto Website](https://quarto.org/) for a more in depth introduction to Quarto docs. The "Document Purpose & Description" section can be deleted when creating a new analysis.## Document Components1. **YAML (yet another markdown language) header** - the code contained at the top of this document bounded by a set of `---` marks. This is a place to store document metadata and instructions for how the rest of the quarto doc should run code chunks and format rendered documents. See here for more information --\>[Authoring Docs](https://quarto.org/docs/authoring/front-matter.html){.uri} and [YAML HTML Options](https://quarto.org/docs/reference/formats/html.html#citation). We currently have the YAML set so when you Render this document a self-contained .html document will be created, this file format and formatting options can be further customized.2. **Code chunks** - code chunks are like mini code scripts. In a quarto doc you can run these individually by clicking on the green arrow/triangle in the upper right corner of the chunk OR using your "run code" shortcut. These chunks are allow to user to write text in between code, and provide instructions on how the code should be incorporated into the rendered document, including how figures should be displayed. See here for more information --\>[Code chunk options (work in both quarto and Rmarkdown docs)](https://quarto.org/docs/computations/execution-options.html){.uri}3. **Markdown** - This bullet list is written in markdown. Markdown is a simple text and formatting language. Use this to write text in the rendered report document. Here is a quick guide to markdown syntax --\>[Markdown cheat sheet](https://www.rstudio.com/wp-content/uploads/2015/02/rmarkdown-cheatsheet.pdf){.uri}. You can also use values and/or objects generated in your code to populate in your markdown text, see here ---\>[How to use a code value in markdown text](https://www.earthdatascience.org/courses/earth-analytics/multispectral-remote-sensing-modis/add-variables-to-rmarkdown-report/){.uri}## Code Chunk Example - Run CodeThe primary function of this document is to run code. Below is an example of a code chunk that runs R code. These are where data cleaning and analysis steps are contained. In Rstudio you can run code several ways:```{r example-code-chunk}#| echo: true# Code chunk example# The coding language is indicated in the very first line, within the curly brackets {}. # Use `#| ` to set options for the individual code chunk below the first line of the code chunk.# A code chunk execute option `echo: true` enables the printing of code output in rendered report # here is example R code. To write comments within a code chunk, use a # at the beginning of the line, lines with # will not be run within a code chunk.test <-1+1test_2 <- test *6test_2```1. line-by-line (just like a normal .R file)2. chunk-by-chunk (press the green right-facing triangle at the top right corner of the code chunk)3. all code top to bottom ()## Render Quarto Docs - Build Report DocThe secondary function of this document is to build a report document. This requires all the code in the Quarto document to run cleanly without errors. To build the report, click the `Render` button shown below. This will start a multistage process that runs all the code and markdown sequentially and formats it based on the YAML header and code chunk execute options. There are two ways to alter/customize the resulting report:1. The YAML header (described above). For example, to set the default execution options for all code chunks, you can add options in the YAML (`execute:`)2. Set run options for individual code chunks. More details here --\><https://quarto.org/docs/computations/execution-options.html>The final template .html product will look something like this:# Analysis TemplateBelow is the code required for eDNA analysis## Setup### Load packagesWe used the `renv` package to manage & record package versions. This records the version of R, package versions, and their dependencies into a file `renv.lock` at the root of this `.Rproj`. The purpose of this is to allow users on different computer setups and in the future to reproduce the results of this project in the same way. This is because different package versions can behave differently version to version, and even introduce errors or produced different results. By using the same exact set of package versions that the authors used to write the code, the same results can be reproduced and prevent unintended errors.`Renv` creates a separate custom project library that is separate from the system library typically used when the `library()` function is called. This means that when you recreate the `renv` project library by calling `renv::restore()` you will likely need to install many packages. Many may be already installed in your system library, but they will need to be installed again into this custom `renv` project library with the exact version recorded in the `renv.lock` file. All files, folders, and scripts within the `.Rproj` will share the `renv` project library.For more information about `renv` and how to add/update the `renv` project library -\>[ ](https://rstudio.github.io/renv/articles/renv.html).::: callout-warningAlthough this package management strategy is meant to reduce errors among future users, it is not perfect. The authors found that the `renv` project library could not yet be recreated on a M1 chip Apple computer (2023-10).:::```{r load-packages}#| output: false# load renv packages from project renv.lock file#renv::restore()library(renv)library(tidyverse) # this include tidyr, ggplot2, dplyr for data manipulation and graphinglibrary(vegan)library(data.table)library(superheat)library(calecopal) # color palettelibrary(viridisLite) # color blind & gray scale friendly color palettelibrary(cowplot) # create graphs with multiple components in a gridlibrary(viridisLite)#library(fantaxtic)library(magrittr) # contains pipe function "%>%"#devtools::install_github("joey711/phyloseq")library(phyloseq)packageVersion("phyloseq")library(owmr)library(devtools)#BiocManager::install("MicrobiotaProcess")library(MicrobiotaProcess)library(patchwork)#devtools::install_github("r-lib/conflicted")library(conflicted)```### Create paths to folders```{r file-directories-paths}# create file paths to folders within the project for easy & consistent reference throughout the analysis# path to raw data folderdir_data_raw <-file.path("./data/raw")# path to processed / clean data folderdir_data_proc <-file.path("./data/processed")# path to code outputs folderdir_products <-file.path("./products")# path to code scripts folderdir_scripts <-file.path("./scripts")# path to images for reportdir_images <-file.path("./images")######## create folders if they do not exist# data folderif (!dir.exists(file.path("./data/"))) {dir.create(file.path("./data/"))}# raw data folderif (!dir.exists(dir_data_raw)) {dir.create(dir_data_raw)}# processed data folderif (!dir.exists(dir_data_proc)) {dir.create(dir_data_proc)}# outputs folderif (!dir.exists(dir_products)) {dir.create(dir_products)}# scripts folderif (!dir.exists(dir_scripts)) {dir.create(dir_scripts)}# images folderif (!dir.exists(dir_images)) {dir.create(dir_images)}```### Set analysis date```{r set-analysis-date}# set this variable here, used throughout analysis for file naming. Used to differentiate between analyses using the same analysis template on different dataanalysis_date <-"2021_11"```### Load data```{r load-data}esv_data_all <-read_csv(file.path(dir_data_raw, "JVB1710-MiFishU-read-data.csv"))# ESV - exact sequence variant# to see a snippet of the data run the following code in the console #glimpse(esv_data_all)```## Clean & save data### Clean & organize dataThe code chunk below "clean-data" is separated into its own `./scripts/clean_data.R` file and sourced by calling `source(file.path(dir_scripts, clean_data.R)` here. If you are looking save space in this analysis document, or want to reuse this exact same data cleaning protocol in different analysis document. The same goes for "create-nearshore-offshore-data" and "create-near-offshore-species-list" chunks, they could be their own scripts or be added to a clean_data script. If you choose to do this, the "example-source-clean-data" chunk below can be used.```{r read-run-clean-data}# only run the script if the output .csv file does not yet exist# load data if it already exsistsif (!file.exists(file.path(dir_data_proc, paste0("esv_clean_data_rv_", analysis_date, ".csv")))) {source(file.path(dir_scripts, "clean_data.R")) # run cleaning script in scripts foldercat("ran cleaning script")} else {esv_clean_data_df <-read_csv(file.path(dir_data_proc, paste0("esv_clean_data_rv_", analysis_date, ".csv")))cat("loaded exsisting clean data")}```### Create Nearshore & Offshore dataframes (or load if they already exist)```{r create-nearshore-offshore-data}# store file names in objects to reuse & reference. Change here if needed, and will be automatically used through out codefile_offshore_esv <-paste0("offshore_esv_", analysis_date, ".csv")file_nearshore_esv <-paste0("nearshore_esv_", analysis_date, ".csv")# if esv file does not exist: run code and save# if esv file does exist: load file##### offshoreif (!dir.exists(file.path(dir_data_proc, file_offshore_esv))) {# create offshore_esv dataframe offshore_esv <- esv_clean_data_df %>%select(2, 6, 7, 8, 9, 10, contains('COS'), contains('AFZ'), contains('DFZ'), contains('NMS'))# save dataframe as .csvwrite_csv(offshore_esv, file =file.path(dir_data_proc, file_offshore_esv))# if file exists - read it in} else (offshore_esv <-read_csv(file.path(dir_data_proc, file_offshore_esv)) )###### nearshoreif (!dir.exists(file.path(dir_data_proc, file_nearshore_esv))) {# create offshore_esv dataframe nearshore_esv <- esv_clean_data_df %>%select(2, 6, 7, 8, 9, 10, contains('COS'), contains('SHR'), contains('LGN'), contains('FRR'))# save dataframe as .csvwrite_csv(nearshore_esv, file =file.path(dir_data_proc, file_nearshore_esv))# if file exists - read it in} else (nearshore_esv <-read_csv(file.path(dir_data_proc, file_nearshore_esv)) )```### Create Nearshore & Offshore species dataframes (or load if they already exist)```{r create-near-offshore-species-list}# store file names for reusefile_offshore_esv_sp <-paste0("offshore_esv_sp_", analysis_date, ".csv")file_nearshore_esv_sp <-paste0("nearshore_esv_sp_", analysis_date, ".csv")# similar to code chunk above - if output file already exists, data will load, if not data created#### offshore speciesif (!dir.exists(file.path(dir_data_proc, file_offshore_esv_sp))) {# create offshore species list offshore_esv_sp <- esv_clean_data_df %>%select(2, 10, contains('COS'), contains('AFZ'), contains('DFZ'), contains('NMS'))# save dataframe as .csvwrite_csv(offshore_esv_sp, file =file.path(dir_data_proc, file_offshore_esv_sp))# if file exists - read it in} else (offshore_esv_sp <-read.csv(file.path(dir_data_proc, file_offshore_esv_sp)))#### nearshore speciesif (!dir.exists(file.path(dir_data_proc, file_nearshore_esv_sp))) {#create nearshore species list nearshore_esv_sp <- esv_clean_data_df %>%select(2, 10, contains('COS'), contains('SHR'), contains('LGN'), contains('FRR'))# save dataframe as .csvwrite_csv(nearshore_esv_sp, file =file.path(dir_data_proc, file_nearshore_esv_sp))# if file exists - read it in} else (nearshore_esv_sp <-read_csv(file.path(dir_data_proc, file_nearshore_esv_sp)))```## Data Analysis### Unique Species```{r unique-species}#| eval: falsen_distinct(esv_clean_data_df$Species)unique(esv_clean_data_df$Species)```### Unique Species - nearshore```{r unique-nearshore-species}#| eval: false# filters out (removes) all rows from the dataframe nearshore_esv that contain only zeros, keeping only the rows that have at least one non-zero elementns_esv <- nearshore_esv[rowSums(nearshore_esv==0, na.rm=TRUE)<ncol(nearshore_esv), ]#filters out all rows from nearshore_esv whose sum across all columns is 0 or negative, keeping only rows where this sum is positive#nearshore_esv[rowSums(nearshore_esv[])>0,]n_distinct(nearshore_esv$Species)unique(nearshore_esv$Species)```### Unique Species - offshore```{r unique-offshore-species}#| eval: falsen_distinct(offshore_esv$Species)unique(offshore_esv$Species)```### Unique Genera```{r unique-genera}#| eval: falsen_distinct(esv_clean_data_df$Genus)unique(esv_clean_data_df$Genus)```### Unique Families```{r unique-families}#| eval: falsen_distinct(esv_clean_data_df$Family)unique(esv_clean_data_df$Family)```### Transform data to Phyloseq Object```{r create-phyloseq-object}#| output: false# AM - Can we find a more reproducible way of selecting columns? Explore what column names are we want to keepnames(esv_clean_data_df)#create an OTU/ESV count tableesv_read_count_m <- esv_clean_data_df %>%select(-c(1,3:20)) %>%# remove first column, and columns 3-20data.frame(row.names =1) %>%# convert result to data.frame object with row names#remove_prefix(c("X")) %>% # remove prefix from column names if needed - remove first # in this line# names() %>% # if prefix removed above - remove first # in this line (could make this into a conditional statement)as.matrix() # convert data.frame to matrix object classESV <- phyloseq::otu_table(esv_read_count_m, taxa_are_rows =TRUE)#create a taxonomy tableesv_read_tax_m <- esv_clean_data_df %>%select(-c(1,3)) %>%# remove columns 1 and 3 - AM will this be standard across datasets? can we make this into positive selection to make more concise?select(-c(9:195)) %>%# remove columns 9 through 195data.frame(row.names =1) %>%# convert to data.frame class with row namesas.matrix() # convert to matrix objectTAX <- phyloseq::tax_table(esv_read_tax_m) #create a sample metadata samp_meta <-read_csv(file.path(dir_data_raw, "SampleMetaData.csv")) # Need to move this file into raw foldersamp_meta_df <-data.frame(samp_meta, row.names =1) # use 1st column as row namessamp_metadata <- phyloseq::sample_data(samp_meta_df) # create phyloseq metadata#create phyloseq objectphyseq <-phyloseq(ESV, TAX)# merge phyloseq metadata with phyloseq objectphyseq1 <-merge_phyloseq(physeq, samp_metadata)```### Explore Components of Phyloseq Object (if needed)```{r explore-phyloseq}#sample_names(physeq1)#rank_names(physeq1)#sample_variables(physeq1)#rank_names(physeq1)```### Subset Phyloseq Samples - Keep Environmental Samples Only```{r phyloseq-env-samples_only}physeq_env <-subset_samples(physeq1, Control_or_Envtl =="Envtl")```### Subset Phyloseq Environmental Samples - ESVs with marine fish and mammal annotations```{r create-env-marine-fish-mammal, echo=FALSE}#Select for only Classes: Actinopteri and Mammalia, which will remove Classes: Aves and Amphibia, as well as the unannotated ESVs physeq_env_mar <- physeq_env %>%subset_taxa(Class %in%c("Actinopteri", "Mammalia")) %>%# select only boney fish and mammalssubset_taxa(!(Order %in%c("Primates"))) %>%# remove primatessubset_taxa(!(Family %in%c("Bovidae", "Canidae", "Felindae", "Hominidae", "Suidae"))) # remove non-marine mammals# Explore contamination in samples# physeq_mar <- subset_taxa(physeq_env, Class %in% c("Actinopteri", "Mammalia"))# physeq_mar# physeq_contam_primate <- subset_taxa(physeq_env, !(Order %in% c("Primates")))# physeq_contam_primate# physeq_contam_husbandry <- subset_taxa(physeq_env, !(Family %in% c("Bovidae", "Canidae", "Felindae", "Hominidae", "Suidae")))# physeq_contam_husbandry```### Environmental Marine Samples With More Than 10,000 Sequences```{r env-marine-more-than-10k}# create new phyloseq object - Environmental Marine fish and mammal samples with less than 10,000 sequencesphyseq_env_mar_lt10k <-prune_samples(sample_sums(physeq_env_mar)<10000, physeq_env_mar)# create new phyloseq object - Environmental Marine fish and mammal samples with greater than 10,000 sequencesphyseq_env_mar_pruned <-prune_samples(sample_sums(physeq_env_mar)>=10000, physeq_env_mar)```### Standardize number of reads to the median sequencing depth```{r std-reads-meadian-seq-depth}total <-median(sample_sums(physeq_env_mar_pruned))standf <-function(x, t=total) round(t *(x /sum(x)))physeq_mar_std <-transform_sample_counts(physeq = physeq_env_mar_pruned, fun = standf)``````{r}# AM - can we remove?#rcurve(physeq_mar_std, add_sample_data = TRUE)#p <- ggrare(physeq_mar_std, step = 500, color = "Zone", label = "Sample", se = FALSE)#sample_data(physeq_mar_std)#TopNESVs <- names(sort(taxa_sums(physeq_mar_std), TRUE)[1:30])#physeq_mar_std_30 <- prune_taxa(TopNESVs, physeq_mar_std$Species)#print(physeq_mar_std_30)```### Subset taxa```{r, subset_taxa}#| echo: false# Chat GPT description - x represents the abundance data for a single taxon (like a species) across all samples. The function calculates the coefficient of variation, which is the standard deviation (sd(x)) divided by the mean (mean(x)), of this abundance data. The criterion for keeping a taxon in the dataset is that this coefficient of variation must be greater than 3.0. In simpler terms, it filters out taxa whose abundance varies less across samples compared to those whose abundance varies more.sf =filter_taxa(physeq_mar_std, function(x) sd(x)/mean(x) >3.0, TRUE)# Extract the taxonomic data as a data frametax_data <-as.data.frame(phyloseq::tax_table(sf))##### Subset by boney fish# Check if 'Actinopteri' is in the 'Class' columnif("Actinopteri"%in% tax_data$Class) {# Subset if 'Actinopteri' is found physeq_std_Act <-subset_taxa(sf, Class =="Actinopteri")cat("Subset created successfully\n")} else {# Handle the case where 'Actinopteri' is not foundcat("'Actinopteri' not found in the Class column\n")}##### Subset by Mammals # Check if 'Mammalia' is in the 'Class' columnif("Mammalia"%in% tax_data$Class) {# Subset if 'Mammalia' is found physeq_std_Mar <-subset_taxa(sf, Class =="Mammalia")cat("Subset created successfully\n")} else {# Handle the case where 'Mammalia' is not foundcat("'Mammalia' not found in the Class column\n")}##### subset by Sharks & Rays# Check if 'Chondrichthyes' is in the 'Class' columnif("Chondrichthyes"%in% tax_data$Class) {# Subset if 'Chondrichthyes' is found physeq_std_Con <-subset_taxa(sf, Class =="Chondrichthyes")cat("Subset created successfully\n")} else {# Handle the case where 'Chondrichthyes' is not foundcat("'Chondrichthyes' not found in the Class column\n")}# create dataframe from phyloseq object#physeq_std_Act_df <- psmelt(physeq_std_Act)```# Visualizing### Barplot - Families - Environmental Samples (excluding Controls)```{r figure-barplot-families_enviro}#| label: Fam-Env-samples#| title: Families from Environmental samples#| fig-width: 6#| fig-asp: 0.618#| out-width: "40%"#| fig-align: center#| fig-format: "png"plot_bar(physeq_env, fill ="Family")```### Barplot - Normalized Environmental Samples & Marine ESVs- Family```{r figure-barplot-normal-env-marine-family}plot_bar(physeq_mar_std, fill ="Family")```- Order```{r figure-barplot-normal-env-marine-order}plot_bar(physeq_mar_std, fill ="Order")```### Heatmap - Normalized Environmental Samples & Marine ESVs- Phylum```{r figure-heatmap-normal-env-marine-phylum}plot_heatmap(physeq_mar_std, taxa.label="Phylum")```### Barplot - Alternative formats- Family```{r}plot_bar(physeq_mar_std, x="Zone", fill ="Family")```- Family facet wrap```{r}plot_bar(physeq_mar_std, x="Zone", fill ="Family") +facet_wrap(~Cardinal_direction)``````{r}#plot_net(physeq_mar_std, maxdist=0.4, color="Cardinal_direction", shape="Near_or_Offshore")```### Richness plots```{r}# plot_richness() function produces a ggplot object and can be ammended just like ggplot# source color palettesource(file.path("./scripts/color_palette.R"))Fig_richness <-plot_richness(physeq_mar_std, measures =c("Chao1", "Shannon"), x="Near_or_Offshore", color="Near_or_Offshore",shape="Near_or_Offshore") +theme_light(base_size =12) +scale_color_manual(values = near_offshore_colors) +labs(color ="Region",shape ="Region") +theme(axis.title.x =element_blank())Fig_richnessggsave(plot = Fig_richness, device ="png", filename =file.path("products/figures/Fig_richness_near_offshore_rv.png"),dpi =600,width =6,height =4,units ="in")``````{r}plot_richness(physeq_mar_std, measures=c("Chao1", "Shannon"), x="Near_or_Offshore", color="Zone")```### Ordinations```{r}physeq_mar_std.ord <-ordinate(physeq_mar_std, "PCoA", "bray")# base plot ordination - other iterations are built on this objectpo <-plot_ordination(physeq_mar_std, physeq_mar_std.ord, type ="samples", color ="Zone", shape ="Near_or_Offshore", title ="Nearshore & Offshore Samples")poggsave(po,device ="png",filename =file.path("products/figures/plot_ordination.jpeg"))``````{r, echo=TRUE}Fig_ord_near_off <- po +facet_wrap(~Near_or_Offshore) +scale_color_manual(values = zone_colors,breaks =c("SHR", "LGN", "FRR", "AFZ", "DFZ", "NMS"),guide =guide_legend(override.aes =list(size =3, shape =c(16,16,16,17,17,17)))) +guides(shape =guide_legend(override.aes =list(size =3, shape =c(1,2)))) +theme_light() +labs(shape ="Region") +coord_fixed()Fig_ord_near_offggsave(plot = Fig_ord_near_off, device ="png", filename =file.path("products/figures/Fig_ord_near_offshore_rv.png"),dpi =600,width =7,height =3.5,units ="in")### USE for Prelim Report``````{r}po +facet_wrap(~Cardinal_direction) +scale_color_manual(values = zone_colors,breaks =c("SHR", "LGN", "FRR", "AFZ", "DFZ", "NMS"),guide =guide_legend(override.aes =list(size =3, shape =c(16,16,16,17,17,17)))) +guides(shape =guide_legend(override.aes =list(size =3, shape =c(1,2)))) +theme_light() +labs(shape ="Region") +coord_fixed()``````{r}Fig_ord_near_off_grid <- po +facet_wrap(~factor(Zone, levels =c("SHR", "LGN", "FRR", "AFZ", "DFZ", "NMS"))) +scale_color_manual(values = zone_colors,breaks =c("SHR", "LGN", "FRR", "AFZ", "DFZ", "NMS"),guide =guide_legend(override.aes =list(size =3, shape =c(16,16,16,17,17,17)))) +guides(shape =guide_legend(override.aes =list(size =3, shape =c(1,2)))) +theme_light() +labs(shape ="Region") +coord_fixed()Fig_ord_near_off_gridggsave(plot = Fig_ord_near_off_grid, device ="png", filename =file.path("products/figures/Fig_ord_near_offshore_grid_rv.png"),dpi =600,width =7,height =5,units ="in")#ORDER BY NEARSHORE & OFFSHORE###USE for Prelim Report```### Marine Actinopteri (Bony Fish) bar plot```{r, cache=TRUE}# reorder s4 phyloseq object Sample data by Near_or_Offshore# data_to_rearrange <- physeq_std_Act@sam_data$Sample# data_to_order_by <- physeq_std_Act@sam_data$Near_or_Offshore# ordering <- order(data_to_order_by)# physeq_std_Act@sam_data$Sample <- data_to_rearrange[ordering]fish_color <-viridis(n =90, option ="turbo")Fig_bony_fish <-plot_bar(physeq_std_Act, fill ="Family") +theme_classic() +theme(axis.title.x =element_blank(),axis.text.x =element_blank(),axis.line.x.bottom =element_blank(),axis.ticks.x =element_blank(),legend.position="bottom") +guides(fill =guide_legend(ncol=10)) +scale_y_continuous(expand =c(0,0)) +scale_fill_manual(values = fish_color)# extract legendbony_fish_legend <- cowplot::get_legend(Fig_bony_fish# + # create some space to the left of the legend# theme(legend.box.margin = margin(0, 0, 0, 0)))# place figure without legend on gridFig_bony_fish_leg <- cowplot::plot_grid( Fig_bony_fish +theme(legend.position="none") )# add legend to girdFig_bony_fish_leg2 <- cowplot::plot_grid( Fig_bony_fish_leg, bony_fish_legend,align ="v",nrow =2,ncol =1,rel_heights =c(3,1),scale =1) #Fig_bony_fish_leg2ggsave(plot = Fig_bony_fish_leg2,device ="png",filename =file.path("products/figures/Fig_bony_fish.png"),dpi =600,width =7,height =5,scale =2,units ="in")Fig_bony_fish_leg2```### Marine Mammals```{r}bar_plot_marine_mammals <-plot_bar(physeq_std_Mar, fill ="Species")#jpeg(file="saving_plot2.jpeg")ggsave(bar_plot_marine_mammals,device ="jpeg",filename =file.path("products/figures/Fig_bar_marine_mammals.png"))plot_bar(physeq_std_Act, fill ="Family")dev.off()```### Marine Chondrichthyes (Sharks & Rays)```{r bar_plot_chondrichthyes}# First, check if the 'physeq_std_Con' object exists in the global environmentif(exists("physeq_std_Con")) {# Check if it's a valid phyloseq object and contains species-level dataif (!is.null(physeq_std_Con) &&"Species"%in%rank_names(physeq_std_Con)) {# Plot if it's a valid phyloseq object with species data p <-plot_bar(physeq_std_Con, fill ="Species") # add ggplot graph things here! print(p) } else {cat("The phyloseq object 'physeq_std_Con' is either null or does not contain species-level data.\n") }} else {# If the object does not exist in the global environmentcat("'physeq_std_Con' does not exist.\n")}```### Top10```{r}OTUnames10 =names(sort(taxa_sums(sf), TRUE)[1:10])sf10 =prune_taxa(OTUnames10, sf)sf10near =names(subset(sample_data(sf), Near_or_Offshore=="Nearshore"))sf10off =names(subset(sample_data(sf), Near_or_Offshore=="Offshore"))top10 <-names(sort(taxa_sums(sf), decreasing =TRUE)[1:10])phyloseq::tax_table(sf)[top10,]print(sf10near)print(sf10off)```###Subsetting and pruning Nearshore and Offshore sample sets```{r}sample_variables(physeq_mar_std)phyloseq::tax_table(sf10)[top10,]sfNearshore <-subset_samples(sf, Near_or_Offshore =="Nearshore")sfNearshoresfOffshore <-subset_samples(sf, Near_or_Offshore =="Offshore")sfOffshore#test <- phyloseq::prune_samples(sample_sums(sfNearshore)>=1, sfNearshore)#test#sfNearshore#test2 <- phyloseq::prune_samples(sample_sums(sfOffshore)>=1, sfOffshore)#test2#sfOffshoresfNearshore <- phyloseq::prune_taxa(taxa_sums(sfNearshore)>=1, sfNearshore)sfNearshoresfOffshore <- phyloseq::prune_taxa(taxa_sums(sfOffshore)>=1, sfOffshore)sfOffshoresfNearshoresfOffshore```###Subset by Zones```{r}sfSHR <-subset_samples(sf, Zone =="SHR")sfSHRsfLGN <-subset_samples(sf, Zone =="LGN")sfLGNsfFRR <-subset_samples(sf, Zone =="FRR")sfFRRsfAFZ <-subset_samples(sf, Zone =="AFZ")sfAFZsfDFZ <-subset_samples(sf, Zone =="DFZ")sfDFZsfNMS <-subset_samples(sf, Zone =="NMS")sfNMSsfSHR <- phyloseq::prune_taxa(taxa_sums(sfSHR)>=1, sfSHR)sfSHRsfLGN <- phyloseq::prune_taxa(taxa_sums(sfLGN)>=1, sfLGN)sfLGNsfFRR <- phyloseq::prune_taxa(taxa_sums(sfFRR)>=1, sfFRR)sfFRRsfAFZ <- phyloseq::prune_taxa(taxa_sums(sfAFZ)>=1, sfAFZ)sfAFZsfDFZ <- phyloseq::prune_taxa(taxa_sums(sfDFZ)>=1, sfDFZ)sfDFZsfNMS <- phyloseq::prune_taxa(taxa_sums(sfNMS)>=1, sfNMS)sfNMS```###Family counts of physeq_mar_std object (sf)```{r}get_taxa_unique(sf, "Family")#93get_taxa_unique(sfNearshore, "Family")#79get_taxa_unique(sfOffshore, "Family")#53```### Genus counts of physeq_mar_std object (sf)```{r}get_taxa_unique(sf, "Genus")#260get_taxa_unique(sfNearshore, "Genus")#232get_taxa_unique(sfOffshore, "Genus")#83```### Species counts of physeq_mar_std object (sf)```{r}get_taxa_unique(sf, "Species")#389get_taxa_unique(sfNearshore, "Species")#346get_taxa_unique(sfOffshore, "Species")#81```### Top```{r}top5sf_species <-sort(tapply(taxa_sums(sf), phyloseq::tax_table(sf) [, "Family"], sum), decreasing =TRUE)[1:5]names(sort(taxa_sums(sf), decreasing =TRUE)[1:10])```### Top Species```{r}devtools::install_github("gmteunisse/fantaxtic")require("fantaxtic")top10sf_species <-top_taxa(sf, tax_level ="Species",n_taxa =10,grouping ="Near_or_Offshore")top10sf_species```### Top Nearshore Species```{r}top10sfNear_species <-top_taxa(sfNearshore, tax_level ="Species",n_taxa =10)top10sfNear_species```### Top Offshore Species```{r}top10sfOff_species <-top_taxa(sfOffshore, tax_level ="Species",n_taxa =10)top10sfOff_species```###To convert phyloseq object to dataframe use metagMisc```{r}#| eval: falsedevtools::install_github("vmikk/metagMisc")library(vmikk/metagMisc)```###Species per Zone```{r}sfSHR_species <-top_taxa(sfSHR, tax_level ="Species",n_taxa =30)sfSHR_species#psmelt(sfSHR_species)#write_csv(sfSHR_species, file.path(dir_data_proc, "sfSHR_species.csv"))``````{r}sfLGN_species <-top_taxa(sfLGN, tax_level ="Species",n_taxa =30)sfLGN_species``````{r}sfFRR_species <-top_taxa(sfFRR, tax_level ="Species",n_taxa =30)sfFRR_species``````{r}sfAFZ_species <-top_taxa(sfAFZ, tax_level ="Species",n_taxa =30)sfAFZ_species``````{r}sfDFZ_species <-top_taxa(sfDFZ, tax_level ="Species",n_taxa =30)sfDFZ_species``````{r}sfNMS_species <-top_taxa(sfNMS, tax_level ="Species",n_taxa =30)sfNMS_species```
There are two ways to alter/customize the resulting report:
